function ret=irf_nad(c1,c2,E,G,r1,r2,sig_x,sig_y,res_x,res_y,coint_x,coint_y,trend_uk,trend_us,m,n,T,p,h,LE,lr,rest,st,SIG1,SIG2,bayes,quadr,kilian,bias_term_F,bias_term_C,tot_process)
% G_x0=G(1:end,1:m);
% m=num of exogenous variables
% n=num of endogenous variables
% p=num of lags;
% n=num of endogenous variables
% trend=1 if domestic variables have a trend in coint, 0 if not
% lr=1if we go by long run restriction;
for i=1:m:m*(p-1)-(m-1)
e{i}=[E(1:end,i:i+m-1)];
end
% e=coefficient matrices for exogenous variables
for i=1:m:m*(p-1)-(m-1)+m
g{i}=[G(1:end,i:i+m-1)];
end
% g=coefficient matrices for endogenous variables-exogenous block
G_en=G(1:end,m*(p-1)+m+1:end);
for i=1:n:n*(p-1)-(n-1)
g_{i}=[G_en(1:end,i:i+n-1)];
end
% g_=coefficient matrices for endogenous variables-endogenous block
ret.g=g;
ret.g_=g_;
ret.e=e;
sig_yx=cell2mat(ret.g(1))*sig_x;
res_reduced_y=sig_yx/(sig_x)*res_x+res_y;
c_y=c2+sig_yx/(sig_x)*c1;
c_tot=[c1;c_y];
% ret.c_tot=c_tot;
ret.res_red_y=res_reduced_y;
total_res=[res_x;res_reduced_y];
ret.total_res=total_res;
% if LE>0
% SIGMA=total_res*total_res'*(1/(T-p+1));
% else
% SIGMA=total_res*total_res'*(1/(T-p));
% end
if LE==1&& m==length(total_res(:,1))-1;
SIGMA=total_res*total_res'*(1/(T-(p-1)-(p-1)*m-1));
ret.SIGMA=SIGMA;
end
if LE==1&& m<length(total_res(:,1))-1;
SIGMA=total_res*total_res'*1/(T-(p-1));
ret.SIGMA=SIGMA;
end
if LE==0&& m==length(total_res(:,1))-1&&r1==0&&r2==0;
    SIGMA=total_res*total_res'*(1/(T-p-(p-1)*m-1));
ret.SIGMA=SIGMA;
end
if LE==0&& m<length(total_res(:,1))-1;
    SIGMA=total_res*total_res'*1/(T-p);
ret.SIGMA=SIGMA;
end
if LE==0&& m==length(total_res(:,1))-1&&r1>0;
    SIGMA=total_res*total_res'*(1/(T-p));
ret.SIGMA=SIGMA;
end
if LE==1&& m==length(total_res(:,1))-1&&trend_us==1;
SIGMA=total_res*total_res'*(1/(T-(p-1)-(p-1)*m-2));
ret.SIGMA=SIGMA;
end
if LE==1&& m==length(total_res(:,1))-1&&trend_us==1&&quadr==1;
SIGMA=total_res*total_res'*(1/(T-(p-1)-(p-1)*m-3));
ret.SIGMA=SIGMA;
end
if bayes==1&& m==length(total_res(:,1))-1;
       SIGMA=[[SIG1;zeros(1,m)] zeros(m+1,1)];SIGMA(m+1,m+1)=1;
% SIGMA=blkdiag(SIG1,SIG2);
% SIG_yx=cell2mat(ret.g(1))*SIG1;
% SIGMA(m+1:end,1:m)=SIG_yx;
% SIGMA(1:m,m+1:end)=SIG_yx';
% SIGMA(m+1:end,m+1:end)=SIG2+SIG_yx*SIG1^-1*SIG_yx';
% SIGMA=SIG1;
          ret.SIGMA=SIGMA;
end
if (bayes==1&&m <length(total_res(:,1))-1)||(bayes==1&&tot_process==1);
%       SIGMA=[[SIG1;zeros(1,m)] zeros(m+1,1)];SIGMA(m+1,m+1)=1;
SIGMA=blkdiag(SIG1,SIG2);
SIG_yx=cell2mat(ret.g(1))*SIG1;
SIGMA(m+1:end,1:m)=SIG_yx;
SIGMA(1:m,m+1:end)=SIG_yx';
SIGMA(m+1:end,m+1:end)=SIG2+SIG_yx*SIG1^-1*SIG_yx';
          ret.SIGMA=SIGMA;
end
% if bayes==1&& m==length(total_res(:,1))-1&&tot_process==0;
%        SIGMA=[[SIG1;zeros(1,m)] zeros(m+1,1)];SIGMA(m+1,m+1)=1;
% % SIGMA=blkdiag(SIG1,SIG2);
% % SIG_yx=cell2mat(ret.g(1))*SIG1;
% % SIGMA(m+1:end,1:m)=SIG_yx;
% % SIGMA(1:m,m+1:end)=SIG_yx';
% % SIGMA(m+1:end,m+1:end)=SIG2+SIG_yx*SIG1^-1*SIG_yx';
% % SIGMA=SIG1;
%           ret.SIGMA=SIGMA;
% end
% if LE<1&&(r1>0||r2>0);
%     SIGMA=total_res*total_res'*1/T;
% ret.SIGMA=SIGMA;
% end
% P=chol(SIGMA);
% if exo>=1
%    P=diag(diag(P));
% else
%     P=P;
% end

t_x=coint_x(:,1);
 t_y=coint_y(:,1);
 t_red_y=t_y+sig_yx/(sig_x)*t_x;
 ret.t_x=t_x;
 ret.t_red_y=t_red_y;
B_1=G(:,m+1:m+m*(p-1))+sig_yx/(sig_x)*E;
B_2=G(:,m+m*(p-1)+1:end);
zero=zeros(m+n);
if trend_us<1
    coint_x=coint_x(:,:);
else
    coint_x=coint_x(:,2:end);
end
if trend_uk<1
    coint_y=coint_y(:,:);
else
    coint_y=coint_y(:,2:end);
end
% if r1<1
%     extended_coint_x=[];
% else
    extended_coint_x=[coint_x zero(1:m,1:n)];
if r2<1
  coint_y=zeros(n,length(extended_coint_x));
end
% end
if r2<1;
    coint_reduced_y=[];
end
if r2>=1&&r1>=1
coint_reduced_y=coint_y+sig_yx/(sig_x)*extended_coint_x;
end
if r2>=1&&r1<1
   coint_reduced_y=coint_y;
end
ret.B_1=B_1;
ret.B_2=B_2;
ret.coint_reduced_y=coint_reduced_y;
ret.extended_coint_x=extended_coint_x;
A_1=E;
if r1<1&&r2>=1
    COINT=[zero(1:m,1:m+n);coint_reduced_y];
end
if r1>=1&&r2>=1
COINT=[extended_coint_x;coint_reduced_y];
end
if r2<1&&r1>=1
  COINT=[extended_coint_x;zero(1:n,1:m+n)];  
end
if r1<1 && r2<1
    COINT=[zero];  
end
% COINT=total coint matrix
ret.COINT=COINT;
for i=1:m:m*(p-1)-(m-1)
e{i}=[E(1:end,i:i+m-1)];
end
for i=1:m:m*(p-1)-(m-1)
g{i}=ret.B_1(1:end,i:i+m-1);
end
for i=1:n:n*(p-1)-(n-1)
g_{i}=ret.B_2(1:end,i:i+n-1);
end
zero_2=zeros(m+n);
e(cellfun(@isempty,e))=[];
for i=1:p-1
d{i}=[cell2mat(e(i)) zero_2(1:m,1:n)];
end
ret.g_reduced=g;
ret.g__reduced=g_;
ret.e_reduced=e;
% for i=1:m:m*(p-1)-(m-1)
% d{i}=[e(i) zero_2(1:m,1:n)];
% end
g_(cellfun(@isempty,g_))=[];
g(cellfun(@isempty,g))=[];
d(cellfun(@isempty,d))=[];
for i=1:p-1
f{i}=[g{i} g_{i}];
end
for i=1:p-1
D{i}=[d{i};f{i}];
end
% D=each lagged coefficent matrix as in jmulti
ret.D=D;
% if r1<1&r2<1
%     L1=eye(m+n)+cell2mat(ret.D(1));
% else
% L1=eye(m+n)+ret.COINT+cell2mat(ret.D(1));
% end
L1=eye(m+n)+COINT+cell2mat(ret.D(1));
for i=2:p-1
L{i}=cell2mat(ret.D(i))-cell2mat(ret.D(i-1));
end
Lp=-cell2mat(ret.D(p-1));
for i=2:p-1
F{i}=cell2mat(L(i));
end
F{1}=L1;
F{p}=Lp;
zero_3=zeros((n+m)*p-(n+m));
 F_=[eye((n+m)*p-(n+m))];
 F_1=[F_ zero_3(:,1:n+m)];
FF=[cell2mat(F);F_1];
R_=cell(1,h);
% DERIVATION OF COMPANION MATRIX (DD) FOR RELEVANT VAR;
% if LE==1;
% for i=1:p-1
% D{i}=cell2mat(ret.D(i));
% end
% zero_4=zeros((n+m)*(p-1)-(n+m));
%  D_=[eye((n+m)*(p-1)-(n+m))];
% else
%     for i=1:p
% D{i}=cell2mat(ret.D(i));
%     end
%     zero_4=zeros((n+m)*p-(n+m));
%  D_=[eye((n+m)*p-(n+m))];
% end
%  D_1=[D_ zero_4(:,1:n+m)];
% DD=[cell2mat(D);D_1];
% ret.COMPANION=DD;
% R_=cell(1,h);
% KILIAN STATIONARITY CORRECTION;
if kilian==1;
    if max(abs(eig(FF)))>1
  bias_corr_F=FF;
bias_corr_C=c_tot;  
else
% FF=mean(FF,3);bias_term_F=FF-COMP_MAT;
bias_corr_F=FF-bias_term_F;
% c_tot=mean(c_tot,3);bias_term_C=c_tot-CONS;
bias_corr_C=c_tot-bias_term_C;
if LE==1;
%      bb=abs(bias_corr_F^10000-bias_corr_F^9999);
while max(abs(eig(bias_corr_F)))>1
   for i=1:100
       x(1,:)=1;x(i+1,:)=x(i,:)-0.01;
       bias_term_F1(:,:,1)=bias_term_F;
              bias_term_C1(:,:,1)=bias_term_C; 
      bias_term_F1(:,:,i+1)=x(i+1,:)*bias_term_F1(:,:,i);
            bias_term_C1(:,:,i+1)=x(i+1,:)*bias_term_C1(:,:,i); 
      bias_corr_F1(:,:,1)= bias_corr_F;
            bias_corr_C1(:,:,1)= bias_corr_C; 
    bias_corr_F1(:,:,i+1)=FF-bias_term_F1(:,:,i+1);
            bias_corr_C1(:,:,i+1)=c_tot-bias_term_C1(:,:,i+1);        
%  b(:,:,i)=abs(bias_corr_F1(:,:,i+1)^10000-bias_corr_F1(:,:,i+1)^9999);
if  max(abs(eig(bias_corr_F1(:,:,i+1))))<=1;
    bias_corr_F=bias_corr_F1(:,:,i+1);
    bias_corr_C=bias_corr_C1(:,:,i+1);
    bias_term_F=bias_term_F1(:,:,i+1);
    bias_term_C=bias_term_C1(:,:,i+1);
break
end
   end
end
else
%    while abs( max(bias_corr_F^10000))>0.0000001
%     bb= abs(bias_corr_F^10000-bias_corr_F^9999);
  while max(abs(eig(bias_corr_F)))>1
   for i=1:100
       x(1,:)=1;x(i+1,:)=x(i,:)-0.01;
       bias_term_F1(:,:,1)=bias_term_F;
              bias_term_C1(:,:,1)=bias_term_C; 
      bias_term_F1(:,:,i+1)=x(i+1,:)*bias_term_F1(:,:,i);
            bias_term_C1(:,:,i+1)=x(i+1,:)*bias_term_C1(:,:,i); 
      bias_corr_F1(:,:,1)= bias_corr_F;
            bias_corr_C1(:,:,1)= bias_corr_C; 
    bias_corr_F1(:,:,i+1)=FF-bias_term_F1(:,:,i+1);
            bias_corr_C1(:,:,i+1)=c_tot-bias_term_C1(:,:,i+1);        
% if abs(max(bias_corr_F1(:,:,i+1)^10000))<=0.0000001;
%  b(:,:,i)=abs(bias_corr_F1(:,:,i+1)^10000-bias_corr_F1(:,:,i+1)^9999);
if  max(abs(eig(bias_corr_F1(:,:,i+1))))<=1;
    bias_corr_F=bias_corr_F1(:,:,i+1);
    bias_corr_C=bias_corr_C1(:,:,i+1);
    bias_term_F=bias_term_F1(:,:,i+1);
    bias_term_C=bias_term_C1(:,:,i+1);
break
end
   end 
   end
end
    end
     FF=bias_corr_F;
c_tot=bias_corr_C;  
end

for i=1:h
R_{i}=FF^i;
end
IRF=cell(1,h);
for i=1:h
IRF{i}=R_{i}(1:n+m,1:n+m);
end
ret.FF=FF;
ret.c_tot=c_tot;
ret.IRF=IRF;
% derivation of impulse response function
% if lr==0&&bayes==0
% P=chol(SIGMA,'lower');
% end
% if bayes==1
%  if bayes==1&& m==length(total_res(:,1))-1;
%      P=chol(SIGMA,'lower');
% %  P=[[P;zeros(1,m)] zeros(m+1,1)];
%  else

    [P,flag]=chol(SIGMA,'lower'); 
    if flag>0
        P=chol(nearestSPD(SIGMA));
    end
%  end
if lr==1
    P=(eye(m+n,m+n)-(sum(cat(4,D{:}),4)))*(chol((eye(m+n,m+n)/(eye(m+n,m+n)-(sum(cat(4,D{:}),4)))*SIGMA*(eye(m+n,m+n)/(eye(m+n,m+n)-(sum(cat(4,D{:}),4))))'),'lower'));
end
for i=1:h
S_IRF{i}=ret.IRF{i}*P;
end
S_IRF=[P S_IRF];
ret.S_IRF=S_IRF;
% derivation of structural shocks
SHOCKS=P\total_res;
 shocks=SHOCKS';
shocks_us=shocks(:,1:m);
% derivation of us shocks
ret.shocks=shocks;
ret.shocks_us=shocks_us;IR_US=cell(1,h+1);IR_UK=cell(1,h+1);
 for j=1:m
for i=1:h+1
IR_US{i}=[ret.S_IRF{i}(j:j,1:m)];
end
IRF_US{j}=reshape(cell2mat(IR_US'),h+1,m);
 end
 for j=1:n
for i=1:h+1
IR_UK{i}=[ret.S_IRF{i}(m+j:m+j,1:m+n)];
end
IRF_UK{j}=reshape(cell2mat(IR_UK'),h+1,m+n);
 end
 for i=1:m
IRF_US1{i}=[IRF_US{1,i}(1,:);diff(IRF_US{1,i})];
 end
for i=1:n
IRF_UK1{i}=[IRF_UK{1,i}(1,:);diff(IRF_UK{1,i})];
end
S_IRF_=cell(1,h-1);
    S_IRF_{1}=S_IRF{1};
for i=2:h
S_IRF_{i}=S_IRF{i}-S_IRF{i-1};
end
 


% IRF_UK and US,shocks are used only for historical decomposition code
% shocks_us are used only for U.S HD WHEREAS shocks for UK
% IRF_US/UK{J} GIVES THE EFFECT OF THE OTHER SHOCKS ON VARIABLE J FOR H PERIODS
% ret.IRF_US/UK CAN BE USED FOR IMPULSE RESPONSE PLOTS
if LE==1
    ret.IRF_US=IRF_US1;
ret.IRF_UK=IRF_UK1;
ret.S_IRF=S_IRF_;
else
ret.IRF_US=IRF_US;
ret.IRF_UK=IRF_UK;
ret.S_IRF=S_IRF;
end
if rest==1;
    S_IRF1=cell(1,h-1);
              S_IRF1{1}=S_IRF{1};
    for i=2:h
   S_IRF1{i}=[S_IRF{i}(1:end-st-1,:);S_IRF{i}(end-st:end,:)-S_IRF{i-1}(end-st:end,:)]; 
    end
   ret.S_IRF=S_IRF1;
end


    
ret.p=P;
d1=diag(P);
non_nor_shocks=diag(d1)*shocks';
ret.non_nor_shocks=non_nor_shocks';
% FOR VD ONLY LINE 136 IS NEEDED, FOR HD LINES 136-142 ARE NEEDED
% HISTORICAL DECOMPOSITION INPUTS
% j and i below gives the variable which the historical decomposition is computed for
% US
%  ret2=irf(E,G,2,0,sig_x,sig_y,res_x,res_y,coint_x,coint_y,0,0,7,1,373,6,400);
% s_us=ret2.shocks_us;
% j=1;
% c_us=ret2.IRF_US{j};
% US-FIRST SAMPLE
%  ret2=irf(E,G,sig_x,sig_y,res_x,res_y,coint_x,coint_y,0,0,6,1,335,6,400);
% s_us=ret2.shocks_us;
% j=1;
% c_us=ret2.IRF_US{j};
% UK
%  ret2=irf(E,G,0,4,sig_x,sig_y,res_x,res_y,coint_x,coint_y,0,0,4,6,370,6,400);
% s_us=ret2.shocks_us;
% s_open=ret2.shocks;
% j=1;
% c_us=ret2.IRF_US{j};
% i=1;
% c_open=ret.IRF_UK{i};
% CANADA
%  ret_=irf(E,G,sig_x,sig_y,res_x,res_y,coint_x,coint_y,0,0,4,6,228,6,250);
% s_us=ret.shocks_us;
% s_open=ret.shocks;
% j=1;
% c_us=ret.IRF_US{j};
% i=1;
% c_open=ret.IRF_UK{i};
% SWEDEN
%  ret_=irf(E,G,sig_x,sig_y,res_x,res_y,coint_x,coint_y,0,0,4,6,159,6,250);
% s_us=ret.shocks_us;
% s_open=ret.shocks;
% j=1;
% c_us=ret.IRF_US{j};
% i=1;
% c_open=ret.IRF_UK{i};